package scu.maqiang.fes;

import static scu.maqiang.numeric.Constants.GP3;
import static scu.maqiang.numeric.Constants.GW3;

import java.util.Arrays;
import java.util.HashSet;
import java.util.stream.IntStream;

import scu.maqiang.mesh.Mesh;
import scu.maqiang.numeric.CG;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

public class FES3H20CTE extends FES{
	public FES3H20CTE(Mesh mesh, FES3H201 fes1, FES3H202 fes2) {
		super(mesh, 4);
		this.fes3H201 = fes1;
		this.fes3H202 = fes2;
	}
	
	/**
	 * @brief 组装不考虑振动与阻尼的直接耦合法的热力耦合整体迭代矩阵
	 * @param constCoef 材料系数向量, 结构均匀各向同性，向量至少有9个分量, 分别表示
	 *                  constCoef[0] - 材料密度rho
	 *                  constCoef[1] - 比热容c
	 *                  constCoef[2] - 热传导系数kxx
	 *                  constCoef[3] - 参考温度 TRef
	 *                  constCoef[4] - 热传导时间theta-迭代格式中的theta值
	 *                  constCoef[5] - 弹性模量 E
	 *                  constCoef[6] - 泊松比 Nu
	 *                  constCoef[7] - 热模量 Beta
	 *                  constCoef[8] - dt 迭代时间步长 
	 * @param tpElasticStiff 组装的弹性刚度矩阵类型
	 * @param tpHeatStiff 组装的热刚度矩阵类型
	 * @param tpHeatMass 组装的热质量矩阵类型
	 * @param tpCoupled 组装的热力耦合矩阵类型
	 * @param LHS 迭代左端矩阵
	 * @param RHS 迭代右端矩阵
	 */
	public void assembleDierctCouplingMatrix_Type1(double[] constCoef, BVPType tpElasticStiff, BVPType tpHeatStiff, 
			                                 BVPType tpHeatMass, BVPType tpCoupled, SRMatrix LHS, SRMatrix RHS) {
		double rho = constCoef[0];
		double c = constCoef[1];
		double kxx = constCoef[2];
		double TRef = constCoef[3];
		double theta = constCoef[4];
		double E = constCoef[5];
		double Nu = constCoef[6];
		double Beta = constCoef[7];
		double dt = constCoef[8];
		double coefM = rho * c / TRef;
		double coefK = kxx / TRef;
		
		int dofHeat = mesh.getnPerElement();
		int dofElastic = mesh.getnPerElement() * mesh.getnDim();
		double[][] ENu = new double[2][dofHeat];
		Arrays.setAll(ENu[0], i -> E);
		Arrays.setAll(ENu[1], i -> Nu);
		double[][] vecBeta = new double[1][dofHeat];
		Arrays.setAll(vecBeta[0], i -> Beta);
		
		double[][] matK = new double[1][dofHeat];
		Arrays.setAll(matK[0], i -> coefK);
		double[][] matM = new double[1][dofHeat];
		Arrays.setAll(matM[0], i -> coefM);
		int nv = mesh.getNv();
		
		double[][] HeatStiff = new double[dofHeat][dofHeat];
		double[][] HeatMass = new double[dofHeat][dofHeat];
		double[][] ElasticStiff = new double[dofElastic][dofElastic];
		double[][] HeatToElastic = new double[dofElastic][dofHeat];
		double[][] ElasticToHeat = new double[dofHeat][dofElastic];
		double[][] eleLMatrix = new double[dofHeat][dofHeat];
		double[][] eleRMatrix = new double[dofHeat][dofHeat];
		int[] idxHeat = new int[dofHeat];
		for(int i = 0, n = mesh.getNt(); i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			int[] idxElastic = getDof(ele, nDim);
			Arrays.setAll(idxHeat, ii -> ele[ii] + nDim * nv);
			FES3H202.elementElasticityStiff(ele, coordRegion, ENu, tpElasticStiff, ElasticStiff);
//			System.out.println(MVO.toString(ElasticStiff));
			elementDirectCouplingMatrix1(ele, coordRegion, vecBeta, tpCoupled, HeatToElastic);
//			System.out.println(MVO.toString(HeatToElastic));
			for(int si = 0; si < dofElastic; si++) {
				for(int sj = 0; sj < dofHeat; sj++) {
					ElasticToHeat[sj][si] = HeatToElastic[si][sj];
					HeatToElastic[si][sj] = -HeatToElastic[si][sj];
				}
			}
			LHS.assemble(idxElastic, ElasticStiff);
			LHS.assemble(idxElastic, idxHeat, HeatToElastic);
			LHS.assemble(idxHeat, idxElastic, ElasticToHeat);
			RHS.assemble(idxHeat, idxElastic, ElasticToHeat);
			
			FES3H201.elementStiff(ele, coordRegion, matK, tpHeatStiff, HeatStiff);
			FES3H201.elementMass(ele, coordRegion, matM, tpHeatMass, HeatMass);
//			System.out.println(MVO.toString(HeatStiff));
//			System.out.println(MVO.toString(HeatMass));
//			System.exit(0);
			for(int si = 0; si < dofHeat; si++) {
				for(int sj = 0; sj < dofHeat; sj++) {
					eleLMatrix[si][sj] = HeatMass[si][sj] + theta * dt * HeatStiff[si][sj];
					eleRMatrix[si][sj] = HeatMass[si][sj] - (1 - theta) * dt * HeatStiff[si][sj];
				}
			}
			LHS.assemble(idxHeat, eleLMatrix);
			RHS.assemble(idxHeat, eleRMatrix);
		}
	}
	
	/**
	 * @brief 组装考虑振动与而不考虑阻尼的直接耦合法的热力耦合整体迭代矩阵
	 * @param constCoef 材料系数向量, 结构均匀各向同性，向量至少有9个分量, 分别表示
	 *                  constCoef[0] - 材料密度rho
	 *                  constCoef[1] - 比热容c
	 *                  constCoef[2] - 热传导系数kxx
	 *                  constCoef[3] - 参考温度 TRef
	 *                  constCoef[4] - 热传导时间theta-迭代格式中的theta值
	 *                  constCoef[5] - 弹性模量 E
	 *                  constCoef[6] - 泊松比 Nu
	 *                  constCoef[7] - 热模量 Beta
	 *                  constCoef[8] - dt 迭代时间步长
	 *                  constCoef[9] - alpha Newmark迭代格式参数
	 * @param tpElasticStiff 组装的弹性刚度矩阵类型
	 * @param tpElasticMass 组装的弹性质量矩阵类型
	 * @param tpHeatStiff 组装的热刚度矩阵类型
	 * @param tpHeatMass 组装的热质量矩阵类型
	 * @param tpCoupled 组装的热力耦合矩阵类型
	 * @param LHS 迭代左端矩阵
	 * @param RHS 迭代右端矩阵
	 */
	public void assembleDierctCouplingMatrix_Type2(double[] constCoef, BVPType tpElasticStiff, BVPType tpElasticMass, 
			BVPType tpHeatStiff, BVPType tpHeatMass, BVPType tpCoupled, SRMatrix LHS, SRMatrix RHS) {
		double rho = constCoef[0];
		double c = constCoef[1];
		double kxx = constCoef[2];
		double TRef = constCoef[3];
		double theta = constCoef[4];
		double E = constCoef[5];
		double Nu = constCoef[6];
		double Beta = constCoef[7];
		double dt = constCoef[8];
		double alpha = constCoef[9];

		double coefM = rho * c / TRef;
		double coefK = kxx / TRef;
		
		int dofHeat = mesh.getnPerElement();
		int dofElastic = mesh.getnPerElement() * mesh.getnDim();
		double[][] matP = new double[1][dofHeat];
		Arrays.setAll(matP[0], i -> rho);
		
		double[][] ENu = new double[2][dofHeat];
		Arrays.setAll(ENu[0], i -> E);
		Arrays.setAll(ENu[1], i -> Nu);
		double[][] vecBeta = new double[1][dofHeat];
		Arrays.setAll(vecBeta[0], i -> Beta);
		
		double[][] matK = new double[1][dofHeat];
		Arrays.setAll(matK[0], i -> coefK);
		double[][] matM = new double[1][dofHeat];
		Arrays.setAll(matM[0], i -> coefM);
		int nv = mesh.getNv();
		double[][] HeatStiff = new double[dofHeat][dofHeat];
		double[][] HeatMass = new double[dofHeat][dofHeat];
		double[][] ElasticStiff = new double[dofElastic][dofElastic];
		double[][] ElasticMass = new double[dofElastic][dofElastic];
		double[][] HeatToElastic = new double[dofElastic][dofHeat];
		double[][] ElasticToHeat = new double[dofHeat][dofElastic];
		double[][] eleLMatrix = new double[dofHeat][dofHeat];
		double[][] eleRMatrix = new double[dofHeat][dofHeat];
		double[][] eleLMatrix6 = new double[dofElastic][dofElastic];
		double[][] eleRMatrix6 = new double[dofElastic][dofElastic];
		int[] idxHeat = new int[dofHeat];
		for(int i = 0, n = mesh.getNt(); i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			int[] idxElastic = getDof(ele, nDim);
			Arrays.setAll(idxHeat, ii -> ele[ii] + nDim * nv);
			FES3H202.elementElasticityStiff(ele, coordRegion, ENu, tpElasticStiff, ElasticStiff);
			FES3H202.elementElasticityMass(ele, coordRegion, matP, tpElasticMass, ElasticMass);
			elementDirectCouplingMatrix1(ele, coordRegion, vecBeta, tpCoupled, HeatToElastic);
			FES3H201.elementStiff(ele, coordRegion, matK, tpHeatStiff, HeatStiff);
			FES3H201.elementMass(ele, coordRegion, matM, tpHeatMass, HeatMass);
			
			for(int si = 0; si < dofElastic; si++) {
				for(int sj = 0; sj < dofHeat; sj++) {
					ElasticToHeat[sj][si] = HeatToElastic[si][sj];
					HeatToElastic[si][sj] = -HeatToElastic[si][sj];
				}
			}
			
			for(int si = 0; si < dofHeat; si++) {
				for(int sj = 0; sj < dofHeat; sj++) {
					eleLMatrix[si][sj] = HeatMass[si][sj] + theta * dt * HeatStiff[si][sj];
					eleRMatrix[si][sj] = HeatMass[si][sj] - (1 - theta) * dt * HeatStiff[si][sj];
				}
			}
			
			for(int si = 0; si < dofElastic; si++) {
				for(int sj = 0; sj < dofElastic; sj++) {
					eleRMatrix6[si][sj] = 1.0 / (alpha * dt * dt) * ElasticMass[si][sj];
					eleLMatrix6[si][sj] = ElasticStiff[si][sj] + eleRMatrix6[si][sj];
				}
			}
			
			LHS.assemble(idxElastic, eleLMatrix6);
			LHS.assemble(idxElastic, idxHeat, HeatToElastic);
			LHS.assemble(idxHeat, idxElastic, ElasticToHeat);
			LHS.assemble(idxHeat, eleLMatrix);
			
			RHS.assemble(idxElastic, eleRMatrix6);
			RHS.assemble(idxHeat, idxElastic, ElasticToHeat);
			RHS.assemble(idxHeat, eleRMatrix);
		}
		
	}
	
	public void assembleHeatToElastic(double beta, BVPType tpCoupled, SRMatrix HTEMatrix) {
		int nv = mesh.getNv();
		int dofHeat = mesh.getnPerElement();
		int dofElastic = mesh.getnPerElement() * mesh.getnDim();
		double[][] vecBeta = new double[1][dofHeat];
		Arrays.setAll(vecBeta[0], i -> beta);
		double[][] HeatToElastic = new double[dofElastic][dofHeat];
		int[] idxHeat = new int[dofHeat];
		for(int i = 0, n = mesh.getNt(); i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			int[] idxElastic = getDof(ele, nDim);
			Arrays.setAll(idxHeat, ii -> ele[ii] + nDim * nv);
			elementDirectCouplingMatrix1(ele, coordRegion, vecBeta, tpCoupled, HeatToElastic);
			HTEMatrix.assemble(idxElastic, ele, HeatToElastic);
		}		
	}
	/**
	   *  计算单元热力耦合矩阵
	 * @param coord 单元坐标数组
	 * @param beta 热模量系数
	 * @param tp 问题类型
	 * @param K 得到的单元热力耦合矩阵
	 */
	public void elementDirectCouplingMatrix1(int[] ele, double[][] coord, double[][] beta, BVPType tp, double[][] K) {
		double[] x = MVO.col(coord, 0);
		double[] y = MVO.col(coord, 1);
		double[] z = MVO.col(coord, 2);
		double[] dx = new double[20];
		double[] dy = new double[20];
		double[] dz = new double[20];
		double[] v = new double[20];
		double xi, et, zt, detJ, wxarea;
		MVO.fill(K, 0.0);
		switch (tp) {
		case COMMON:
			for(int i = 0; i < 3; i++) {
				xi = GP3[i];
				for(int j = 0; j < 3; j++) {
					et = GP3[j];
					for(int k = 0; k < 3; k++) {
						zt = GP3[k];
						detJ = CG.ISOMap3D(CG::H20ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
						double coefBeta = MVO.dot_product(beta[0], v);
						wxarea = GW3[i] * GW3[j] * GW3[k] * detJ;
						for(int ki = 0; ki < 20; ki++) {
							for(int kj = 0; kj < 20; kj++) {
								K[3 * ki][kj] += coefBeta * wxarea * dx[ki] * v[kj];
								K[3 * ki + 1][kj] += coefBeta * wxarea * dy[ki] * v[kj];
								K[3 * ki + 2][kj] += coefBeta * wxarea * dz[ki] * v[kj];
							}
						}
					}
					
				}
			}
			break;
		case CONSTITUTE_MATRIX_COMMON:
			for(int i = 0; i < 3; i++) {
				xi = GP3[i];
				for(int j = 0; j < 3; j++) {
					et = GP3[j];
					for(int k = 0; k < 3; k++) {
						zt = GP3[k];
						detJ = CG.ISOMap3D(CG::H8ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
						double coefBeta1 = MVO.dot_product(beta[0], v);
						double coefBeta2 = MVO.dot_product(beta[1], v);
						double coefBeta3 = MVO.dot_product(beta[2], v);
						wxarea = GW3[i] * GW3[j] * GW3[k] * detJ;
						for(int ki = 0; ki < 20; ki++) {
							for(int kj = 0; kj < 20; kj++) {
								K[3 * ki][kj] += coefBeta1 * wxarea * dx[ki] * v[kj];
								K[3 * ki + 1][kj] += coefBeta2 * wxarea * dy[ki] * v[kj];
								K[3 * ki + 2][kj] += coefBeta3 * wxarea * dz[ki] * v[kj];
							}
						}
					}
				}
			}
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
	}
	
	/**
	  *  对稀疏矩阵施加温度边界条件
	 * @param L
	 * @param heatBdLabel
	 */
	public void applyBCHeat_MBN(SRMatrix L, int... heatBdLabel) {
		HashSet<Integer> bdSet = mesh.extractBoundaryNodes(heatBdLabel);
		int nv = mesh.getNv();
		bdSet.forEach(label -> L.setElement(label + nDim * nv, label + nDim * nv, 1.0e30));
	}
	
	public void applyBCHeat_MBN(double[] RHS, double val, int... heatBdLabel) {
		HashSet<Integer> bdSet = mesh.extractBoundaryNodes(heatBdLabel);
		int nv = mesh.getNv();
		bdSet.forEach(label -> RHS[label + nDim * nv] = val * 1.0e30);
	}
	
	public void applyBCDisp_MBN(SRMatrix L, Direct direct, int... heatBdLabel) {
		HashSet<Integer> bdSet = mesh.extractBoundaryNodes(heatBdLabel);
		int[] dof = getDof(bdSet, direct, nDim);
		//System.out.println(Arrays.toString(dof));
		IntStream.of(dof).forEach(idx -> L.setElement(idx, idx, 1.0e30));
	}
	
	public void applyBCDisp_MBN(double[] RHS, Direct direct, double val, int... heatBdLabel) {
		HashSet<Integer> bdSet = mesh.extractBoundaryNodes(heatBdLabel);
		int[] dof = getDof(bdSet, direct, nDim);
		IntStream.of(dof).forEach(idx -> RHS[idx] = val * 1.0e30);
	}
	
	public void extractUVH(double[] RHS, double[] temp, double[][] dispUV) {
		for(int i = 0, n = mesh.getNv(); i < n; i++) {
			temp[i] = RHS[i + nDim * n];
			for(int j = 0; j < nDim; j++) {
				dispUV[j][i] = RHS[nDim * i + j];
			}
		}
	}
	
	public FES3H201 fes3H201;
	public FES3H202 fes3H202;
}

